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Gravitational waves emitted during the coalescence of binary neutron star systems are self¬ 
calibrating signals. As such, they can provide a direct measurement of the luminosity distance 
to a source without the need for a cross-calibrated cosmic distance-scale ladder. In general, how¬ 
ever, the corresponding redshift measurement needs to be obtained via electromagnetic observations 
since it is totally degenerate with the total mass of the system. Nevertheless, Fisher matrix studies 
have shown that, if information about the equation of state of the neutron stars is available, it is pos¬ 
sible to extract redshift information from the gravitational wave signal alone. Therefore, measuring 
the cosmological parameters in pure gravitational-wave fashion is possible. Furthermore, the huge 
number of sources potentially observable by the Einstein Telescope has led to speculations that the 
gravitational wave measurement is potentially competitive with traditional methods. The Einstein 
Telescope is a conceptual study for a third generation gravitational wave detector which is designed 
to yield 10® — 10^ detections of binary neutron star systems per year. This study presents the first 
Bayesian investigation of the accuracy with which the cosmological parameters can be measured 
using information coming only from the gravitational wave observations of binary neutron star sys¬ 
tems by Einstein Telescope. We find, by direct simulation of 10® detections of binary neutron stars, 
that, within our simplifying assumptions, Hq, Ulrn, Da, Wq and wi can be measured at the 95% level 
with an accuracy of ~ 8%, 65%, 39%, 80% and 90%, respectively. We also find, by extrapolation, 
that a measurement accuracy comparable with current measurements by Planck is possible if the 
number of gravitational wave events observed is 0(10®“^).We conclude that, while not competitive 
with electro-magnetic missions in terms of significant digits, gravitational wave alone are capable of 
providing a complementary determination of the dynamics of the Universe. 


I. INTRODUCTION 

The family of second generation interferometers Ad¬ 
vanced LIGO [1] began its operations in the last quarter 
of 2015 [2]. Advanced Virgo [3] is scheduled to join the 
LIGO network in 2017, with KAGRA [4] and LIGO 
India [5] to follow afterwards. The detection of grav¬ 
itational waves from the coalescence of merging black 
holes [6-8] led already to important scientific measure¬ 
ments as tests of general relativity [8, 9] and astrophysics 
[8, 10, 11]. Given the expected number of yearly detec¬ 
tions [8, 12, 13], the expectations on the scientific deliv¬ 
erables are high: tests of the strong field dynamics of 
General Relativity [8, 9, 14-16]; a “cosmic distance scale 
ladder”-free determination of the Hubble constant [17- 
19]; a determination of the neutron star equation equa¬ 
tion of state [20-23]. 

Detectors beyond the second generation are already 
being envisaged. For instance, the Einstein gravitational 
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wave Telescope (ET) [24] is a proposed underground de¬ 
tector consisting of three 10 km arm-long Michelson in¬ 
terferometers in a triangular topology with opening an¬ 
gles of 60 degrees [25]. The strain sensitivity is estimated 
as factor 10 better than second generation detectors, 
down to frequencies of 1 — 3 Hz depending on the actual 
configuration of the instrument [26]. The high sensitivity 
promises the detection of a very large number of gravi¬ 
tational waves (GW) signals with large signal-to-noise 
ratios (SNR), thus allowing for unprecedented popula¬ 
tion studies as well as extremely accurate measurements 
of the physical parameters of coalescing binary systems 
[24]. 


A. Cosmological inference with gravitational waves 

When estimating the parameters of GW sources, and 
in particular the coalescences of binary neutron stars and 
black holes, the luminosity distance can be observed di¬ 
rectly [27, 28]. This makes GW an ideal laboratory to 
place samples in the Hubble diagram in a manner that is 
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free from the potential systematics affecting electromag¬ 
netic (EM) methods. Unfortunately, in the vast majority 
of the cases, the redshift cannot be measured from GW 
alone and this piece of information needs to be extracted 
by means other than GW. 

In the recent years, various proposals have been put 
forward to overcome this difficulty. For instance, one 
can assume that the coalescences of compact binaries are 
the progenitors of short gamma ray bursts (sGRB) [29]. 
In this case, coincident observations of a GW event and 
the correspondent sGRB would allow the measurement of 
the luminosity distance from GW and the redshift from 
spectroscopy of the host galaxy[17, 30-32]^. For second 
generation interferometers, this method indicates a rela¬ 
tive accuracy on the measurement of the Hubble constant 
Hq of few percent in the case where 10-15 such events are 
detected. However, whether the coalescences of compact 
binaries are the progenitors of sGRBs is still a matter 
of debate. Also, the fraction of GW events also observ¬ 
able as sGRBs might be as low as 10“^ [35] due to sGRB 
beaming effects. 

An alternative approach, following broadly the argu¬ 
ment first given in Ref. [27], would be to statistically 
identify the possible host galaxies of a GW event to ob¬ 
tain a distribution of possible redshifts associated with 
each GW detection. This method should yield ^ 5% 
percent accuracy on Hq using 20-50 events [18] observed 
by Advanced LIGO/Virgo. A similar methodology has 
also been applied to space-based detectors [36, 37]. 

A few methods aim at extracting the redshift using 
GW observations alone. For example, one can use the 
knowledge of the (rest frame) mass function of NS and 
the measured (redshifted) mass to infer the redshift of 
the source [19, 38]. In this framework, second generation 
interferometers should infer the Hubble constant Hq with 
~ 10% accuracy using about 100 events [19]. 

The results of advanced interferometers can be greatly 
improved by third generation instruments such as ET. In 
fact, ET can probe regions of the Universe where the ef¬ 
fects of the dark energy will be substantial, thus allowing 
an independent sampling of the cosmic history. 

The potentialities of ET have already been investigated 
by various groups [30, 31, 39], concluding that, when only 
a limited set of cosmological parameters is considered, the 
accuracy of the inference is comparable to that of current 
EM measurements. 


B. Outline 

In this paper, we will expand on the approach pro¬ 
posed by Messenger and Read [40] in which if one of the 


^ Kilonovae are also expected EM counterparts to BNS coales¬ 
cences [33]. However their utility as cosmological probes is yet 
unclear due to their intrinsically faint luminosities [e.g. 34] which 
limits the distance at which they can be confidently detected. 


two components is a NS, information about the equa¬ 
tion of state (EOS) allows a direct measurement of the 
rest-frame masses and thus of the source redshift [40]. 
Using Fisher matrix formalism, the authors estimate the 
accuracy with which z can be measured to be ~ 8 — 40%, 
depending on the EOS and on the distance to the source. 
Recently, a similar investigation was carried out in [41] 
using a more realistic Monte Garlo data analysis method. 
The authors concluded that the average uncertainty is 
closer to 40% for a hard EOS and essentially indepen¬ 
dent of redshift. 

Nevertheless, given the large number of sources that 
can be observed by ET and the possibility of combining 
information across them, even the large uncertainty re¬ 
ported in Ref. [41] could be sufficient to obtain interesting 
indications on the accuracy with which ET will measure 
the cosmological parameters. In this paper, we explore 
this idea in a simplihed scenario and conclude that ET 
can indeed set bounds that are comparable to current 
EM measurement. We are interested in the cosmological 
information that can be inferred exclusively from the ob¬ 
servation of gravitational waves. We will thus not discuss 
the potential of coincident GW-EM detections which are 
presented elsewhere [31, 39]. We note here that, because 
of the co-location of the three ET interferometers and 
because of its topology, its expected sky resolution is ex¬ 
tremely poor. Gonsequently, the probability of a success¬ 
ful EM-GW association is a priori very small. Note that 
at the time of ET, second generation detectors are ex¬ 
pected to be operational with improved sensitivities [42]. 
For a substantial fraction of the loudest GW events, the 
sky localisation from a network made of ET and advanced 
detectors will be vastly improved compared to ET alone. 
In this case, some of the aforementioned EM-I-GW meth¬ 
ods might become feasible and used to yield constraints 
on the cosmological parameters. 

The rest of the paper is organised as follows. In foot¬ 
note 1 we cast the problem in a Bayesian framework, and 
identify the necessary components to arrive at the cosmo¬ 
logical inference. In Eq. (10) we describe the procedures 
of simulating GW events and the detector noise, and the 
implementation of the analysis. In Fig. 3 we present the 
results of our simulations and Hnally in table HI we sum¬ 
marise and discuss our results. The mathematical solu¬ 
tion to the problem of the inference of the cosmological 
parameters in the presence of a detection threshold is 
given for completeness in Appendix A. 


II. METHOD 

In this section we present a Bayesian solution to the 
problem of computing posterior probability density func¬ 
tions for a set of cosmological parameters from GW data. 
We broadly follow the presentation in [18]. Note that the 
treatment is not specific to the case of ET. 
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A. Inference of the cosmological parameters in the 
absence of a detection threshold 

Consider a catalogue of GW events £ = {ei,...,e 7 v}- 
Each event is defined as a stretch of data di{t) given by 
the sum of noise ni{t) and a gravitational wave signal 
i.e. 

Ci : di{t) = ni{t) + hi{Qi;t), ( 1 ) 

where 0^ indicates the set of all parameters of the signal 
i. 

The noise is taken to be a stationary Gaussian process 
with a zero mean and covariance defined by its one-sided 
spectral density £'„(/) such that 

p(„.|I)ocexp{-ly^"d/4W^}, 

ot exp{-i (n| n)| (2) 

where I represents all the relevant information for the in¬ 
ference problem, a tilde represents the Fourier transform, 
and where we have introduced a scalar product between 
two real functions A{t) and B{t) as 

, 3 ) 

The likelihood of observing the event is then given by 

p(G|0i,S,I) oc exp|-^ {d, -hi\di- hi)| (4) 

where S is the signal model that relates the signal param¬ 
eters Qi to a gravitational wave signal h. Moreover, the 
signal-to-noise ratio (SNR) p can be succinctly written 
as 

P= Vih\h). (5) 

The posterior distribution for any parameter in our signal 
model S is related to the likelihood in Eq. (4) through the 
application of Bayes’ theorem 

p(0j|ei, S, I) cx _p(0i|S, T)p{ei\Qi, S, I) (6) 

where p(0i|S,I) is the prior probability distribution for 
the parameters 0^. When multiple independent detec¬ 
tors are included in the analysis, the likelihood (Eq. (4)) 
generalises to 

p(e,|0„S,I)=[]p(ef |0„S,I). (7) 

k 

For this work, we are only interested in the pos¬ 
terior probability for a subset of parameters O = 


{Ho,£ljnj£lA, ■ ■ ■}■ Therefore, we marginalise over the 
remaining subset of parameters 9i, i.e. 

S, I) = J d9i p(0i|ej, S, I) 

= J d0ip(0j,f1|S,I)p(ei|0i,O, S,I) 

= p(f2|S, I) IA9,p{9,\a, S, l)p(e,|0l, O, S, I) 

= p(f2|S, !)/:(£„ f2), (8) 

Where we have introduced the so-called “quasi¬ 
likelihood” [43] 

C{ei,n) = yd0jp(0i|f1,S,I)p(ei|0i,O, S,I). (9) 

Finally, the posterior for given an ensemble of events 
£ can be shown to be 

p{^}\£,s,l) = pin\S,l)l[£{e„n). ( 10 ) 

i 

Therefore, in order to obtain the posterior for f2, we need 
to perform a multi-dimensional integral in Eq. (9) for 
each of the GW events. The description of this procedure 
and the generation of data follow in Eq. (10). 

III. ANALYSIS 

In this section we describe the simulation that was per¬ 
formed. Firstly, we outline the generation of the data, 
consisting of the GW signal model, the astrophysical and 
cosmological assumptions regarding the source popula¬ 
tion, and the simulation of the detector noise. Secondly, 
we show the data analysis implementation with which 
the simulated data was analysed. In particular, we de¬ 
scribe the construction of the quasi-likelihood, and it sub¬ 
sequent use to arrive at our cosmological inference. The 
GW signals, and the detector noise have been generated 
using the LIGO Analysis library (LAL) [44]. 

A. Astrophysical and cosmological assumptions 

The NS masses are distributed according to a Gaussian 
distribution with a mean of 1.35M0 and a standard devi¬ 
ation of O.15M0 [45] which is assumed constant through¬ 
out the cosmic history. For the NS equation of state 
we consider three cases; a hard EOS, a medium and a 
soft EOS. They are labelled as MSI [46], H4 [47] and 
SQM3 [48]. We investigate these three cases since in [40] 
it was shown that the accuracy with which the redshift 
can be measured depends on the magnitude of the phys¬ 
ical effects related to the details of the EOS. One can 
think of these three cases as an optimistic, a realistic 
and a pessimistic one, respectively. 
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The events are distributed uniformly in the cosine of 
the inclination, polarisation and time of arrivals. The 
events are also uniformly distributed in comoving vol¬ 
ume. Their redshifts are sampled from the probability 
density given by [49] 


piz\n) 


dR{z) 1 
dz 4?(Zniax) 


( 11 ) 


where R{z) is the cosmic coalescence rate. It is worth 
nothing that p{z\^) is an explicit function of Q. The 
differential cosmic coalescence rate is equal to 


dR{z) dV roe{z) 

dz dz 1 + z 


where rg is the local rate, e(z) is the cosmic star forma¬ 
tion rate and V is the comoving volume. In a Friedmann- 
Robertson-Walker-Lemaitre (FRWL) universe, the rate 
of change of V with z is given by 


dz (l+z)2iJ(z)’ 


(13) 


where we have introduced the Hubble parameter 


JI(z) = + Z)^ + IIfc(I -I- z)2 -I- rtAE{z,w(z)) 

(14) 


and the luminosity distance [50] 

DLiz)=\ (1 + ^)/o"r^ for r!fc =0 

[ for llfc < 0 

(15) 

Hg is the Hubble constant, 17^ is the matter fractional 
density, IIa is the fractional energy density of dark en¬ 
ergy, flfc = 1 — is the curvature. Finally 

£^(Z, W{z)) = (1 -f 2)3(l-H»o-H»i)g-3»iz/(l-Hz) 

is a convenient parametrisation to capture the effects of 
the redshift evolution of dark energy [51]. For (l we chose 
fiducial values of 

(h, fl™, flA, Hfc, wg, u;i)“ = (0.7,0.3,0.7,0,-1,0), (17) 

where h = iJo/lOOkm-Mpc^r.s^r Even though FT hori¬ 
zon distance is ~ 37 Gpc (z ~ 4.15 for our fiducial cos¬ 
mology), we limit our analysis to Zmax = 2 as this corre¬ 
sponds approximately the sky averaged horizon distance 
of 13 Gpc for BNS systems [52]. For simplicity, we de¬ 
cided to assume a star formation rate e(z) that does not 
change with redshift and is therefore irrelevant for our 
problem. 

We simulated 1,000 binary NS events as observed by 
ET. The parameters Oi of each individual source have 
been generated according to the assumptions described in 



FIG. 1. Network SNR distribution of the 1,000 BNS events 
generated sampling Eqs. (11) and (13) for a our fiducial cos¬ 
mology Eq. (17). 


Eq. (10). The corresponding waveform h{f, ©i) was then 
added to Gaussian noise which is coloured according to 
the amplitude spectral density shown in Fig. 2. In Fig. 1 
we show the network SNR distribution in the ET detector 
computed using Eq. (5). 

Differently from most existing literature, we do not hl- 
ter our sources with any SNR threshold. If we were to do 
so, we would be introducing a selection bias [53]. Note 
that, due to the potentially large number of sources ob¬ 
served by ET and their distribution in co-moving volume, 
the vast majority of them will in practice not be detected 
in a search which uses the SNR as decision statistics. It is 
known that ignoring these unregistered sources leads to a 
significant bias in the estimation of “global” parameters, 
see [53]. 

The main reason for the emergence of biases is inti¬ 
mately linked to the functional form for the prior on 
z Eq. (11). Since Eq. (11) quantihes the prior expec¬ 
tation regarding the distribution of sources in the co¬ 
moving volume, it is an explicit function of fl. The quasi¬ 
likelihoods for the majority of our simulated events are 
almost uniform in D, see Section HIG2, therefore our 
inference is greatly influenced by the prior distribution: 
if one were to analyse sources that are louder than some 
threshold SNR the overall population of events would 
appear on average closer than the actual population. At 
the same time, the observed distribution of z would fol¬ 
low the actual cosmological distribution. Since p{z\^) 
Eq. (11) relates Dl and z via D, this leads for estimates 
of Dm —)■ 1 and h —)■ 0. Similarly, if one were to consider 
only events that are quieter than a given SNR thresh¬ 
old, Dm —> 0 and h —>■ 1, thus Da —>-1. In Appendix 
A, we give a mathematical solution to the problem of 
inferring D that accounts for sets of unobserved events. 
However, the solution in A is not computationally treat¬ 
able with current techniques and for a large number of 
events, therefore we opted for an SNR selection threshold 
of 0 and analysed all simulated events. Furthermore, our 
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choice relies on the capacity of distinguishing between 
low SNR GW signals and low SNR background events 
due to noise in the detector. A discussion of this prob¬ 
lem can be found in [53]. The triangular configuration of 
ET provides an additional tool to study the distribution 
of signal and noise low SNR events. Thanks to its topol¬ 
ogy, ET admits the construction of a null stream which is 
devoid of any GW signal as the sum of the outputs of the 
individual Michelson detectors [25]. Being a pure noise 
process, the analysis of the null stream can be used to 
understand the SNR distribution of noise events which 
can then be used to infer the SNR distribution of quiet 
sources. 


B. The signal model 


In the previous paragraph we introduced the signal 
model S without specifying its properties. In this section, 
we lay out the assumptions that go in the construction 
of S. 

In modelling the GW from a binary system, we limit 
our analysis to the inspiral phase of the coalescence pro¬ 
cess. We model the inspiral using an analytical frequency 
domain 3.5 post-Newtonian waveform in which we ig¬ 
nore amplitude corrections and the effects of spins. This 
is not a big limitation as NS are expected to be slow 
rotators[54]. In particular, we use the so-called TaylorF2 
approximant [55], which can be written as 

h(0;/) = A(0)/-"/V^(®-/), (18) 

where the waveform is written in terms of the amplitude 
A(0) and the phase $(0;/). 


The amplitude of the waveform A(0) is given by 


^(0) oc — Qit, Ip, a, S) (19) 

where we have introduced the chirp mass A4 = 
l(m\ + 1712 )^^^, Dl is the luminosity distance 
defined in Eq. (15), {a, 6) signify the sky position of the 
source, and (t, ^p) give the orientation of the binary with 
respect to the line of sight [55]. 

The wave phase can be written in the form 


^(0;/) = 27r/G -(pc- 


In/] /(”-5)/3, 


( 20 ) 


n—0 


where the ipn are the so-called post-Newtonian coeffi¬ 
cients (see e.g. [56]), which are functions of the com¬ 
ponent masses mi and m 2 , and (G, (pc) are the time and 
phase of coalescence. Note that all masses are defined 
in the observer frame, and the rest frame mass m^est is 
related to the observed mass through 


w = mrest(l + 2 ) , (21) 

where z is the redshift of the GW source. 

The description of the phase in Eq. (20) assumes that 
the object is a point particle, and thus cannot be tidally 
deformed. However, since we consider all of our events 
are binary NS coalescences, we modify the gravitational 
wave phase in Eq. (20) by including the Hnite-size contri¬ 
butions to the phase. These in turn depend on the tidal 
deformability A(mrest) [21] of the star which is a func¬ 
tion of its equation of state and its rest frame mass. The 
finite-size contributions to the GW phase, as a function 
of observed masses, is given by 


d>tidai(0; /) = 


a—1 


3Ao(l -l- z)^ 
1287? Afs 




(3179 - 919 xa - 2286 x" + 260 x^) (ttM/)' 


( 22 ) 


where the sum is over the components of the binary, Xo = 
ma/M, Xa = X{ma) where mo are the component masses, 
M is the total mass, and ry = mim 2 /M^. 


Knowledge of the EOS and using information encoded 
in the GW tidal phase contribution allows to measure 
the redshift of the source [40]. While the EOS is not 
known yet, various studies have shown that it could be 
possible to infer it from observations of BNS with second 
generation detectors [20-23, 57, 58]. In what follows, we 
will assume that the nature of the NS interior is known. 


C. Data analysis 

For our analysis, we assumed a noise curve for ET cor¬ 
responding to the “B” configuration [59], corresponding 
to the projected sensitivity achievable with the current 
technologies Fig. 2. Given the anticipated rates of com¬ 
pact binary coalescences [13], the detection rates of bi¬ 
nary NS systems in ET are expected to he in the range 
10® - 10^ yr-®[24]. 

The parameters 9 of our signal model are the compo¬ 
nent masses mi and m 2 , inclination t, polarisation 
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Hz 

FIG. 2. Amplitude spectral density for ET in the “B” config¬ 
uration. 


component masses. For the common parameters we used 
uniform probability distributions on the 2-sphere for sky 
position (a, S) and orientation (ip, l) and uniform in the 
time of coalescence tc with a width of 0.1 seconds around 
the actual coalescence time. For the first marginalisa¬ 
tion, we choose uniform sampling distributions for both 
Dl and z in the intervals [1,10®] Mpc and [0,2], respec¬ 
tively. 

For the component masses, the priors were different 
across the different runs; each EOS in fact predicts not 
only the functional form of the tidal deformability A(m) 
that enters in the phase of the GW waveform, but also 
the maximum permitted mass of the NS itself. Therefore, 
for the three EOSs under consideration we used maxi¬ 
mum expected rest frame mass Mmax of 2.8, 2.0, 2 .OM 0 
for MSI, H4 and SQM3 respectively. The prior prob¬ 
ability distribution for the component masses was then 
chosen to be uniform between IMq and Mmax- 


right ascension a and declination 6 , the time of coales¬ 
cence tc, the phase of coalescence (pc luminosity distance 
Hi and the redshift z. In our analysis we ignored the 
presence of spins, as it is believed to be small in binary 
NS systems [54]. We analysed our ensemble of sources 
assuming that the EOS of the NS is known, thus ac¬ 
counting to a total of three analysis runs (one for each of 
our predefined hard, medium and soft EOSs). To obtain 
the posterior probability distribution on the cosmologi¬ 
cal parameters we proceeded in three steps. Firstly, 
we analysed each source to compute a quasi-likelihood as 
a function of the redshift z and the luminosity distance 
D^. Secondly, these quasi-likelihoods are then converted 
into quasi-likelihoods as a function of the cosmological 
parameters as shown in Eq. (9). Einally, the poste¬ 
rior probability function for the cosmological parameters 
given an ensemble of events are computed from Eq. (10). 


1, Obtaining the quasi-likelihood 
Eor each event Cj, we compute the quantity 

C(ei, Dl, z,9.) = J dXp{X\Q, S, I)p(ei|A, H, S, I) (23) 

that is a partially marginalised quasi-likelihood, where 
the marginalisation is done on all parameters that are 
not relevant to the inference of O. These are A = 
{mi,m2,ip,t,(pc,tc,a,5). The further marginalisation 
over z and Dl will be described later on. For the time be¬ 
ing, let us describe the details of the analysis for the com¬ 
putation of Eq. (23). The above integral was computed 
using a Nested Sampling algorithm [60] implemented sim¬ 
ilarly to what described in [61]. For each of the three 
analysis runs, we chose the same prior probability dis¬ 
tributions for all parameters, with the exception of the 


2. Cosmological inference 


The marginalisation over the redshift and luminosity 
distance was then performed as follows: once a cosmolog¬ 
ical model is introduced, z and Dl are not independent 
parameters anymore, but they are related unequivocally 
by Eq. (15), thus - after some algebra which can be found 
in [18] - we are left with the following integral to compute 


C(ei,n) = 


dzp{z\D,Sl)C{ei,DL(z),z,D), (24) 


where p(z|0. S'!) is given in Eq. (11) and we chose, con¬ 
sistently with the sources generation, Zmax = 2. 

One of the problems we needed to overcome in order 
to perform the integral in Eq. (24) was how to represent 
C{ei, Dl(z), z,rt) in a tractable way. In fact, one of the 
outputs of the Nested Sampling algorithm is a set of sam¬ 
ples drawn from the integrand in Eq. (23) which is dif¬ 
ficult to manipulate - in particular difficult to integrate 
- without making any assumptions about the underlying 
probability distribution. 

A possible treatment of the problem would be to use 
the samples from Eq. (23) and approximate it using a nor¬ 
malised histogram. This procedure was successfully used 
in other unrelated studies [14], however, for our purposes 
it is not accurate enough. In fact, an histogram represen¬ 
tation is dependent on a parameter, the bin size (or equiv¬ 
alently the number of bins once the range is specified), 
which cannot be inferred from the data but has to be 
chosen arbitrarily. The majority of the quasi-likelihoods 
in Eq. (24) tend to be very uniform over the cosmological 
parameter space for individual sources and, as noted in 
Section III A, the inference is strongly dominated by the 
prior on z. Most sources are close or below the detection 
threshold of the detector, thus a single source, in general, 
yields very little information about the underlying cos¬ 
mology, therefore any small fluctuation in the histogram 
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approximation due to the random variation of the num¬ 
ber of samples in any specific bin would be amplified and 
would eventually lead to a biased estimate of the poste¬ 
rior probability density for Q. As an example, compare 
the panels in Fig. 3. The left panel shows samples from 
Eq. (24) for a source having a network SNR=25. In this 
case, the isoprobability contours are almost consistent 
with a normal distribution. The right panel shows in¬ 
stead a source with network SNR=4. In this case the 
samples are almost uniformly distributed, thus if one 
were to approximate Eq. (24) with a two dimensional 
histogram, different choices of the bin size would result 
in different approximations, which would yield to very 
different inference of ft. Instead, we decided to follow a 
different course of action. Given a set of samples, for any 
partition of the parameter space, the resulting probabil¬ 
ity distribution of the observations is always a multino¬ 
mial distribution, therefore the “probability distribution” 
of the occurrences in each bin is a Dirichlet distribution. 
The above property defines a Dirichlet Process [62] which 
can be used to define an analytical representation of the 
underlying probability distribution of which we have only 
a finite number of samples available. For the mathemat¬ 
ical details and definitions, the reader is referred to the 
original paper by Ferguson [62] or to the more recent 
discussions in [63]. We used the approximate variational 
algorithm [64] as implemented in [65] to find the Dirichlet 
Process Gaussian Mixture Model to represent the inte¬ 
grand in Eq. (24). The output of this procedure is an 
analytical representation of the target probability distri¬ 
bution as an infinite mixture of Gaussian distributions 
which is analytical and continuous. This form can then 
be used as the integrand in Eq. (24) and the integral it¬ 
self can be evaluated using another Nested Sampling al¬ 
gorithm whose output can then be combined to compute 
the posteriors for Eq. (10). 

For this last marginalisation, we assume priors on D 
that are uniform in all parameters. In particular, h G 
[0.1,1.0], e [0.0,1.0], Ha € [0.0,1.0], wo € [-2,0] 
and finally mi € [—1,1]- 

IV. RESULTS 

For each of the equations of state assumed 
we considered various scenarios within the set of 
(/i, Dm, Ha, mo, mi). We show posterior distributions 
for the cosmological parameters obtained from the joint 
analysis of 1, 000 events. We also report confidence inter¬ 
vals from number of events greater than 1, 000 obtained 
via extrapolation. 

Moreover, since all EOS yield very similar results, we 
choose to report only posteriors for the cosmological pa¬ 
rameters obtained from the MSI equation of state. 

We computed posterior distribution functions for three 
distinct cases: 

i. flat FRWL universe. Fig. 4; 


ii. general FRWL universe. Fig. 5; 

iii. general FRWL universe, with Dark Energy parame¬ 
ters, Fig. 6. 

In the following subsections we report on the values 
of the various cosmological parameters in each different 
cosmological model we considered. We note here that, 
as expected, the accuracy of the cosmological parame¬ 
ters measurement is better for the flat case and get pro¬ 
gressively worse with the increasing dimensionality of the 
model under consideration. Also, all uncertainties we re¬ 
port are at the 95% confidence level. 

Posterior distributions on the parameters of all cosmo¬ 
logical models from the analysis of 10^ BNS events are 
reported in Figs. 4, 5 and 6. 

Depending on the actual astrophysical rate, ET will 
observe between 10^ and 10^ BNS events per year [24]. 
It is computationally unfeasible at the moment to simu¬ 
late and analyse in a realistic way such a large number 
of events. We therefore extrapolated the results we ob¬ 
tained for 1,000 sources to the maximum expected num¬ 
ber of events. We assumed that the central limit theorem 
holds, in other words that our posteriors for 10^ events 
are approximately Gaussian, and simply scaled the vari¬ 
ance of the one dimensional posteriors by the number of 
sources N. Tables I,II and III show the extrapolation 
of the 95% widths {2a) to a number of events ranging 
from 10^ to 10^ for all relevant parameters for each of 
the cosmological models we considered in our analysis. 

A. Ho 

We find that 10^ BNS observations yield the following 
results: for the model (i) the accuracy is 0.05(7%) which 
remains approximately constant in the general case (ii) 
and worsens to 0.08(11%) in the general case (iii). In 
comparison with other GW studies, we find our mea¬ 
surements to be significantly worse. For instance [18] 
reports a 95% accuracy of ~ 10% with second genera¬ 
tion interferometers and similarly so do [32] and [38]. In 
comparison to traditional methods, the most constrain¬ 
ing measurement from the Planck experiment in conjuc- 
tion with other methods reports an accuracy of ~ 0.5% 
[66]. So with 10^ GW sources current measurements are 
far more accurate than those obtained with our method 
using ET. However, when extrapolating to the potential 
number of observable BNS systems, we find the Planck- 
like accuracy is reached with ^ 10® BNS observations. A 
further order of magnitude improvement is seen for 10^ 
BNS observations, see Table 1. 

B. Dm and Ha 

We find that with 10® BNS events can be mea¬ 
sured with an accuracy of 0.125(47.5%), 0.135(45%) and 
0.19(65.5%) for the models (i), (ii) and (iii) respectively. 





FIG. 3. Sample marginal likelihoods (Eq. (23)) for 2 and Dl- The left panel shows the marginal likelihood for an event with 
a network SNR=25. The right panel shows the marginal likelihood for an event with a network SNR=4. On both panels, the 
lines indicate the source distance and redshift. 
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N 

Model''''^^ 


10® 


lO'^ 


10® 


10® 


10" 

Flat FRW 

0.5 

X 10"® 

1.6 

X 

10“® 

0.5 

X 

10"® 

1.6 

X 

10"® 

0.5 

X 

10"® 

General FRW 

4.6 

X 10"® 

1.5 

X 

10“® 

4.6 

X 

10"® 

1.5 

X 

10"® 

4.6 

X 

10"® 

General FRW-kDE 

0.8 

X 10"® 

2.5 

X 

10“® 

0.8 

X 

10"® 

2.5 

X 

10"® 

0.8 

X 

10“® 


TABLE I. 95% accuracies on the measurement the reduced Hubble parameter for various detected numbers of sources for the 
general 5 parameter case. For 10^ sources, the widths have been computed using our Nested Sampling algorithm, otherwise 
the widths are the result of an extrapolation. 


The above numbers compare very unfavorably with the 
~ 2% yield by Planck[67]. However, if more sources 
are observed the current accuracy is reached with 10®“^ 
sources, depending on the actual cosmological model. 

The situation is similar for the measurement of the cos¬ 
mological constant Ha. With 10® BNS observations we 
find an accuracy of 0.23(37.5%) and 0.275(39%) for the 
models (ii) and (hi) respectively. As a comparison Planck 
reports ~ 0.9% [66]. However, a similar uncertanty is 
reached with 10®“^ sources. The results are summerised 
in Table H. 


C. Dark Energy parameters 

In the case in which the — z relation is modified 
to allow for a time varying cosmological constant, see 
Eq. (16), the parameters wq and wi were included in the 
model. From the analysis of 10® events, at 95% confi¬ 
dence, we find Awq = 0.8 and Awi =0.9. In compari¬ 
son, assuming 10® BNS events with optical counter parts 
and electromagnetic priors on the remaining cosmological 
parameters[31] finds Awq « 0.1 and Awi = 0.3. Most 
recent determination by Planck [67] of the parameters wq 
and Wi reports Awq « 0.2 and Awi « 0.5. Extrapolation 
to the potential number of sources observable in 1 year 


shows that the accuracy on wq and wi can be improved 
by 2 orders of magnitude, thus 1 order of magnitude bet¬ 
ter than the current best estimates, see Table HI. It is 
worth noting that the posterior distributions for wq and 
wi are not very Gaussian, therefore the extrapolations to 
a large number of events might not be as reliable as for 
the other cosmological parameters. 

V. DISCUSSION 

In this study we investigated the potentialities of BNS 
observations with ET as cosmological probes. In particu¬ 
lar, we quantified the cosmological information that can 
be extracted from pure GW observations of BNS. The 
ingredient that allows the measurement of the redshift is 
the knowledge of the NS equation of state and thus of 
the NS tidal deformability. 

We simulated 1,000 events and relied on extrapolation 
to the expected 10^—10^ events per year. The main result 
of this study is that from GW alone, ET could measure 
all cosmological parameters with an accuracy that is com¬ 
parable with current state-of-the-art measurements from 
EM missions. 

This is the very first study of this kind, and therefore it 
should be regarded as a proof-of-principle. Our analysis 
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AD™ 




N 


10® 


lO'^ 

Flat FRW 

1.3 

X 10"" 

4.0 

X 

10“® 

General FRW 

1.3 

X 10"" 

4.2 

X 

10“® 

General FRW-kDE 

1.9 

X 10"" 

0.6 

X 

10“® 


10 ® 

1.3 X 10"^ 
1.3 X 10"^ 
1.9 X 10"^ 


10 ® 

4.0 X 10"® 
4.2 X 10"® 
0.6 X 10"^ 


10 " 

1.3 X 10“® 
1.3 X 10“® 
1.9 X 10“® 


N 

Model'^''^^ 

ADa 

10® 

lO'^ 

10® 

10® 

10" 

General FRW 2.3 x 10"" 

General FRW-I-DE 2.8 x 10“" 

0.7 X 10“" 
0.9 X 10"® 

2.3 X 10"® 
2.8 X 10"® 

0.7 X 10"® 
0.9 X 10"® 

2.3 X 10“® 
2.8 X 10“® 


TABLE II. 95% accuracies on the measurement the matter energy density 11™ and the cosmological constant Qa for various 
detected numbers of sources for the general 5 parameter case. For 10® sources, the widths have been computed using our Nested 
Sampling algorithm, otherwise the widths are the result of an extrapolation. 


N 

ModeT"'^-^ 

Awo 

10® 

lO'^ 

10® 

10® 

10" 

General FRW-I-DE 0.8 x 10® 

2.5 X 10“" 0.8 X 10"" 

2.5 X 10"® 0.8 X 10“® 

N 

Modei^^^^ 

Awi 

10® 

lO'^ 

10® 

10® 

10" 

General FRW-I-DE 0.9 x 10® 

2.9 X 10“" 0.9 X 10"" 

2.9 X 10"® 0.9 X 10“® 


TABLE III. 95% accuracies on the measurement the Dark Energy parameters wo and wi for various detected numbers of 
sources. For 10® sources, the widths have been computed using our Nested Sampling algorithm, otherwise the widths are the 
result of an extrapolation. 


relies on a set of simplifying assumptions which should be 
progressively relaxed for a comprehensive investigation. 
We neglected the NS spins and eventual precession of the 
orbital plane. We do not expect this to be a serious lim¬ 
itation, since NS are expected to be slow rotators [54], 
however it is not clear how a small degree of precession 
would impact the analysis of ET data. There is evidence 
that for binary black holes and NS black hole binaries ac¬ 
counting for precession of the orbital plane leads to more 
accurate distance measurements [68]. If a redshift mea¬ 
surement can be associated to these classes of sources, 
the determination of Q would improve substantially, also 
thanks to the considerably larger volume that is observ¬ 
able by ET. 

Due to its low frequency sensitivity, BNS signals in the 
local Universe {z < 1) will be in the ET sensitive band 
for a time scale of hours to days, depending on the ob¬ 
server frame masses. It is clear that this problem cannot 
currently be investigated using a fully realistic simula¬ 
tion since the generation of the most accurate waveforms 
and cutting edge data analysis algorithms are not yet 
sufficiently fast. A further complication arises from the 
number of signals itself; a detection rate of lO" events 
per year implies an average time delay between signals of 
~ 3 seconds. Given the duration of the signals in band. 


it follows that several signals would be present simulta¬ 
neously in the ET data stream at any given time. No 
systematic study nor algorithm yet exists to investigate 
this problem. 

We further assumed perfect knowledge of the NS EOS. 
While there are indications that second generation inter¬ 
ferometers could measure the EOS [20-23], it is unlikely 
that we would know it without any form of uncertainty. 
However, we do not consider this a serious limitation as 
long as the error bar on any given value of the NS mass 
m and its tidal deformability X(m) is sufficiently small 
to avoid confusion between different EOS. Moreover, we 
did not include in our analysis the merger part of the 
waveform, which, especially for the most distant sources, 
would be in the sensitive band of ET. The inclusion of 
the merger in the analysis would yield more information 
about the BNS mass and spin parameters as well as intro¬ 
ducing a possible further constraint on its redshift [69], 
allowing for more precise measurements. 

A further element that deserves future investigation is 
the effect of detection thresholds. We give a formal solu¬ 
tion to the problem of the inference of the cosmological 
parameters in Appendix A. However, we did not explore 
the details of its practical implementation, which we de¬ 
fer to a future study. 
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-ffQ/100km-s ^ Mpc ^ 


FIG. 4. Posterior distributions for h and Q.m obtained from 
the analysis of 1,000 BNS events for a flat universe (life = 0) 
and no dark energy equation of state. In the one dimensional 
posteriors the dashed lines indicate the 2.5%, 50% and 97.5% 
confidence levels. In the two dimensional posterior distribu¬ 
tion we show the 68%, 95% and 99% confidence regions. On 
all panels the solid (blue) lines indicate the fiducial values. 


We ignored the effects of the detector calibration un¬ 
certainties over the inference of the GW event parameters 
as well their impact over the global inference of O. At the 
end of the first observing run of Advanced LIGO, typical 
amplitude uncertainties (which is relevant for the deter¬ 
mination oi Dl) and phase uncertainties (relevant for the 
estimation of z) were estimated at ~ 10% and 10 degrees, 
respectively [70]. Simulations indicate that ignoring the 
presence of such calibration errors does not lead to signifi¬ 
cant bias in the estimation of the GW parameters, as long 
as the SNR is not very large [71]. However, data analysis 
models for GW analysis that marginalise over calibration 
uncertainties are now available [72] and routinely utilised 
for the actual analysis [10]. The additional calibration 
uncertainty increase the statistical uncertainty of the in¬ 
ferred parameters by a similar amount. In our case, the 
largest source of uncertainty would come from the ampli¬ 
tude calibration and thus on the determination of for 
the GW sources. Assuming, naively, an ET uncertainty 
budget of ~ 10%, we estimate a similar degradation of 
our inference over H. 

We also ignored the effects of weak lensing. Weak lens- 
ing is a zero mean process [73], thus, when averaging over 
thousands of sources, it will not induce an overall bias in 
the estimate of H. A proper account of the lensing uncer¬ 
tainty, would lead to similar conclusions as the detector 
calibration uncertainty. 


Even with the caveats discussed above, our study 
shows that even considering exclusively information com¬ 
ing from GW alone (with no input from any EM asso¬ 
ciation) ET is capable to fully probe the evolution of 
the Universe and determine the value of U with reason¬ 
able accuracy. We emphasise that our results apply to 
a pure GW-based inference of U. A more accurate de¬ 
termination of U from GW-EM joint detections may be 
possible, thus the results presented in this study should 
be regarded as a lower-limit to what the actual poten¬ 
tiality of ET is as a cosmological probe. Nevertheless, we 
showed that GW alone can be a feasible complementary 
and cross-validating route to probe the dynamics of the 
Universe. 
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Appendix A: Inference of the cosmological 
parameters in presence of censored data 

We want to infer the value of the cosmological param¬ 
eters U = (iJo, • ■ ■) given a set of gravitational 

waves observations. Consider a catalogue of gravitational 
wave events E = {ei,..., cat}. Each event is defined as a 
stretch of data di {t) given by the sum of noise rn {t) and a 
gravitational wave signal /ii(0;f), where 0 indicates the 
set of all parameters of the signal and such that the SNR 
p is greater than a given threshold pth- 

The likelihood to observe the event is given by 
p(eij0, S, Pth) I) where S is the signal model that relates 
the signal parameters 0 to a gravitational wave signal 
h. The posterior distribution for the parameters in our 
signal model S comes from the application of Bayes’ the¬ 
orem 

p(0|ei,S,pth,I) (xp(01S,I)p(eil0,S,pth,I) (Al) 

where p(0jS,I) is the prior probability distribution for 
the parameters 0. Given a certain cosmic coalescence 
rate R{(l, z), there will be a certain number M = M{(l, z) 
of gravitational wave events that will not be registered as 
events since their SNR will be below the selected thresh¬ 
old. Nevertheless, they encode information regarding the 
Universe and they must be taken into account in our in¬ 
ference. The likelihood for a non-detected event 
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FIG. 5. Same as Fig. 4, but for a general FRWL universe. 


is given by 

(Efc :Ah) =pie^\n,Rin,z),pth,l) (A2) 

Pth 

(A3) 

For a set of M non-detected events, the likelihood will 
be given by 

M 

(e", Pth) = n ’ ^th) (A4) 

k^l 

= (Cfe ,pth)]^ (AS) 



The number of non-detected events M(il,z) is a nui¬ 
sance parameters which depends on f2, the rate i?(r2, z), 
the detection threshold pthj the observation time T and 
the observed volume pth) as 


M(0,z) =i?(fl,z)t/(i2,pth)T-lV. (A6) 


We are now in the position of writing the posterior dis¬ 
tribution for the cosmological parameters 17: 
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FIG. 6. Same as Fig. 4, but for a general FRWL universe + DE parameters. 


N, Pth, s, I) oc p{^\S, I) 




N 


^ [C^. {e^ , Pth)] piM\n,R{n,z),pth) ■ 

(A7) 


i=l M=0 


It is interesting to verify that Eq. (A7) reduces to also 
Eq. (10) in the limit of pth 0. In this limit we have 


M ^ 0 

p{M\(},R{n,z),pth) -t Sm,o 


(A8) 

(A9) 


p{ef,, pk\i^, R{n, z), pth,l)dpk -> 0 (AIO) 
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therefore the term 


['^fc(«fc>Ah)] 1 


Assuming that the rate z) is given by the integral 
of Eq. (12), we recover the form of the likelihood (10) 
(All) which we used in our study. 


and the non-detection part of the likelihood reduces to 


i^M.o = 1 • (A12) 

M=0 
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